#!/usr/bin/perl -w

use strict;
use autodie;

use Getopt::Std;

use File::Spec;


my $tools_path = '/mnt/ilustre/app/medical/tools';
my $data_path = '/mnt/ilustre/app/medical/.data';


my $var = ". " . $tools_path . "/.var";
system($var);


my %opts =('i', "brca44");

getopts('i:f:nc:', \%opts);

my $ana = ".ana";
my $log = ".log";
my $error = ".error";
my $option = ".opt";
my $chk = ".chk";
my $dir = $ARGV[0];
my $chk_with_parent = File::Spec->catfile($dir, $chk);
open LOG, ">>$log";
open STDERR, ">>$error";
select LOG;

my($intv, $cutadp, $cnvcal, $cnvcal_ref, %sap_name, @sapl_name, %lane_num, %lane_series, $com, $cnvcal_sample);
my $reads_list0 = "reads.list";
my $time_rand = `date +%s%N`;
my $reads_list = $time_rand . ".list";
my $reads_list_new;
my $brca12 = $data_path . "/intervals/1/brca.qiagen.NGHS-102X-Covered.b37.bed";
my $brca44 = $data_path . "/intervals/1/brca.qiagen.NGHS-001X-Covered.b37.bed";
my $cancer50 = $data_path . "/intervals/1/a_cancer_genes_50.bed";

if ( defined $opts{f}) {
	if ($opts{f} =~ /\.bed$/) {
		$intv = $opts{f};
	}else {
		# $intv = undef;
		die "Error: invalid interval file, should be a bed file.\n"
	}

}elsif ($opts{i} eq "brca12") {
	$intv = "$brca12";
	
}elsif ($opts{i} eq "brca44") {
	$intv = "$brca44";
}elsif ($opts{i} eq "cancer50") {
	$intv = "$cancer50";
}else {
	# $intv = undef;
	die "Error: invalid option -i input.\n";
}

# if ($opts{n} == 0) {
	# $cutadp = 0;
# }elsif ($opts{n} == 1) {
	# $cutadp = 1;
# }else {
	# die "Error: invalid option -a input.\n";
# }

if ($opts{n}) {
	$cutadp = 0;
}else {
	$cutadp = 1;
}

if (defined $opts{'c'}) {
	if ($opts{'c'} =~ /.bam/) {
		$cnvcal = 1;
		$cnvcal_ref = $opts{'c'};
		
		$cnvcal_ref = File::Spec->rel2abs( $cnvcal_ref );
		
		
		print LOG "cnv calling ref: $cnvcal_ref\n";
		
	}else {
		die "Error: the reference data of a cnv calling must be a bam file.\n";
	}

}else {
	$cnvcal = 0;
}


print LOG "interval: $intv\n";
print LOG "cut adaptor: $cutadp\n";
print LOG "cnv calling: $cnvcal\n";

# my @ke = keys %opts;
# my @va =values %opts;
# print "@ke\n";
# print "@va\n";
# print "$opts{'c'}\n";



my $abs_path;

if (File::Spec->file_name_is_absolute( $dir )) {

	$abs_path = $dir;
}else {
	$abs_path = File::Spec->rel2abs( $dir );
}


open READS, ">$reads_list";
opendir my $dirh, $dir or die "Error: cannot open $dir: $!";


foreach (readdir $dirh) {
	# my $i = 0;
	if (/^([a-zA-Z0-9\-]+)_[ATGCatgc]+_(L[\d][\d][\d])_R1_([\d][\d][\d])\.(fastq|fq)(\.gz)?/ ||
	/^([a-zA-Z0-9\-]+)_[sS]\d+_(L[\d][\d][\d])_R1_([\d][\d][\d])\.(fastq|fq)(\.gz)?/) {
		my ($nam, $lane, $series) = ($1, $2, $3);
		$sap_name{$1}++;
		$lane_num{$2}++;
		$lane_series{"$2s$3"}++;
		# if (/$lane/) {
			
			# $lane_series[$i] = $3;
		# }
		
		# $i++;
		# my ($lane_num0, $read_num);
		# if ($2 =~ /([1-9]+)/ {
			# $lane_num0 = $1;
		# }else {
			# die "Error: irregular reads name.\n"
		# }
		
		# if ($3 =~ /([1-9]+)/ {
			# $lane_num0 = $1;
		# }else {
			# die "Error: irregular reads name.\n"
		# }
		
		# ($read_num = $3) =~ /[1-9]+/;
		# ($lane_num0 = $2) =~ /[1-9]+/;
		
		my $file = File::Spec->catfile(($abs_path), $_);
		
		print LOG "\$3: $3\n";
		print READS "$file\t";
		s/_R1/_R2/;
		$file = File::Spec->catfile(($abs_path), $_);
		
		print READS "$file\t$nam\t${lane}s$series\n";	
	}
}



closedir $dirh;

close READS;

@sapl_name = keys %sap_name;

if (!-e $chk_with_parent) {

	unlink "$reads_list" and die "Error: The reads directory must has one and only one sample's reads data. the name of reads must be style like NAME_CAGATC_L001_R1_001.fastq[.gz].\n"
		if @sapl_name != 1;
}





my @nam = keys %lane_num;
my @num = values %lane_num;
my @series = keys %lane_series;
print LOG "@nam\n";
print LOG "@num\n";





my $script_path = $tools_path . '/script/';

# my $pre = $script_path . 'pre.sh';
my $samtools= $tools_path . '/samtools-1.2/samtools';
my $cutadaptor = $script_path . 'cutadaptor0702.sh';
my $align = $script_path . 'align0703.sh';
my $realign = $script_path . 'realign.sh';
# my $merge_part = $script_path . 'merge_part.sh';
# my $merge_lane = $script_path . 'merge_lane.sh';
my $recal = $script_path . 'recal0702.sh';
my $cal = $script_path . 'cal.sh';
my $ann = $script_path . 'anno_0424.sh';
my $summ = $script_path . 'summ.sh';
my $fm_top = $script_path . 'format_top_0508.sh';
my $fm_jxl = $script_path . 'format_jxl.sh';
my $quandico = $script_path . 'quandico_0529.sh';



# system $pre;

# -s exist and not zero
if (-s $chk_with_parent) {

	close LOG;
	close STDERR;
	# print LOG "$dir\n";
	chdir $dir;
	
}else {
	my ($sec, $min, $hour, $day, $mon, $year, @last) = localtime;
	$mon +=1;
	$year +=1900;

	my $dir_name = $sapl_name[0] . "_h" . $hour . "m" . $min . "s" . $sec . "_m" . $mon . "d" . $day . "_" . $year;
	$dir_name = lc $dir_name;
	print LOG "$dir_name\n";
	mkdir $dir_name, 0755 or die "cannot make $dir_name directory: $!";

	close LOG;
	close STDERR;
	rename $log, "${dir_name}/$log";
	rename $reads_list, "${dir_name}/$reads_list";
	# rename $ana, "${dir_name}/$ana";
	rename $option, "${dir_name}/$option";
	chdir $dir_name; 

}

$reads_list_new = $sapl_name[0] . "_" . $reads_list0;

rename $reads_list, $reads_list_new; # if many samples run in the same directory and if one sample's run state is Q, then the second run will override the first reads_list if the name is equal.

$reads_list = $reads_list_new;

open LOG, ">>$log";
open STDERR, ">>$error";
select LOG;
open CHK, ">>$chk";





	
open READS1, "$reads_list";

while (<READS1>) {
	chomp;
	my @reads = split /\t/;
	my $reads1 = $reads[0];
	my $reads2 = $reads[1];
	my $sap_name = $reads[2];
	my $base = $reads[3];
	print LOG "\n";
	print LOG "$reads1\n";
	print LOG "$reads2\n";
	print LOG "$base\n";
	print LOG "\n";
	
	# ali sam->bam sort index realign
	# system 'ali.sh $sap_name $reads1 $reads2 $base';
	
	if ($cutadp) {
		print LOG "cut adaptor...\n";
		my @cutadaptor = ($cutadaptor, '-p', $base, $reads1, $reads2);
		# system $cutadaptor, '-p', $base, $reads1, $reads2;
		(system @cutadaptor) ==0 or die "error: @cutadaptor, $?";
		print CHK "${base}.1.cut.fastq\t${base}.2.cut.fastq\n";
		
		print LOG "align...\n";
		my ($reads1, $reads2) = ("$base" . ".1.cut.fastq", "$base" . ".2.cut.fastq");
		my @align = ($align, '-s', $sap_name, '-p', $base, $reads1, $reads2);
		# system $align, '-s', $sap_name, '-p', $base, $reads1, $reads2;
		(system @align) ==0 or die "error: @align, $?";
		print CHK "${base}.sam\n";
		print CHK "${base}.sort.bam\n";
		
		print LOG "realign...\n";
		my @realign = ($realign, '-p', $base);
		(system @realign) ==0 or die "error: @realign, $?";
		print CHK "${base}.realn.bam\n";
	}else {
		
		my @align = ($align, '-s', $sap_name, '-p', $base, $reads1, $reads2);
		# system $align, '-s', $sap_name, '-p', $base, $reads1, $reads2;
		(system @align) ==0 or die "error: @align, $?";
		print CHK "${base}.sam\n";
		print CHK "${base}.sort.bam\n";
		my @realign = ($realign, '-p', $base);
		(system @realign) ==0 or die "error: @realign, $?";
		print CHK "${base}.realn.bam\n";
	
	}
}

close READS1;



#merge by lane
while (my ($lane, $num) = each %lane_num) {
		
	my $com;
	my $realn_suffix = ".realn.bam";
	my $recal_suffix = ".realn.recal.bam";
	my $realn_index_suffix = ".realn.bai";
	my $merge_name = $lane . $realn_suffix;
	my $merge_index = $lane . $realn_index_suffix;
	my $recal_name = $lane . $recal_suffix;
	print LOG "merge name: $merge_name\n";
	
	if ($num > 1) {
		
		print LOG "merge part...\n";
		
		my @samtools_merge = ($samtools, "merge", "-f", $merge_name);
		my @samtools_index = ($samtools, "index", $merge_name);
		
		foreach (@series) {
			if (/$lane/) {
				$com .= ("$_" . $realn_suffix . " ");
				push @samtools_merge, ("$_" . "$realn_suffix");
			}
		}
		$com =~ s/ $//;
		print LOG "merge part inputs: $com.\n";
		print LOG "@samtools_merge.\n";
		
		(system @samtools_merge) ==0 or die "error: @samtools_merge, $?";
		(system @samtools_index) ==0 or die "error: @samtools_index, $?";
		rename "$merge_name.bai", $merge_index;
		# system $samtools, 'merge', '-f', $merge_name, $com;
		# system $samtools, 'index', $merge_name, $merge_index;
		print CHK "${merge_name}\n";
		my @recal = ($recal, '-p', $lane);
		(system @recal) ==0 or die "error: @recal, $?";
		print CHK "${recal_name}\n";
		
		
	}else {
		
		print LOG "no need to merge part.\n";
		my $realn_index;
		foreach (@series) {
			if (/$lane/) {
				$com = "$_" . $realn_suffix;
				$realn_index = "$_" . $realn_index_suffix;
				last;
			
			}
		}
		print LOG "$com.\n";
		rename $com, $merge_name;
		rename  $realn_index, $merge_index;
		print CHK "${merge_name}\n";
		my @recal = ($recal, '-p', $lane);
		(system @recal) ==0 or die "error: @recal, $?";
		print CHK "${recal_name}\n";
	}
}
	
	

	
#merge by sample
my $suffix_name = ".realn.recal.bam";
my $recal_index_suffix = ".realn.recal.bai";
my $merge_name = "1" . $suffix_name;
$cnvcal_sample = $merge_name;
my $merge_index = "1" . $recal_index_suffix;

if (@nam > 1) {

	my $com;
	my @samtools_merge = ($samtools, "merge", "-f", $merge_name);
	my @samtools_index = ($samtools, "index", $merge_name);
	foreach (@nam) {
		$com .= ("$_" . $suffix_name . " ");
		push @samtools_merge, ("$_" . "$suffix_name");
	}
	
	$com =~ s/ $//;
	# push @merge, $com;
	print LOG "merge lane inputs: $com.\n";
	print LOG "merge lane output: $merge_name\n";
	
	(system @samtools_merge) ==0 or die "error: @samtools_merge, $?";
	(system @samtools_index) ==0 or die "error: @samtools_index, $?";
	rename "$merge_name.bai", $merge_index;
	# system $samtools, 'merge', '-f', $merge_name, $com;
	# system $samtools, 'index', $merge_name, $merge_index;
	print CHK "${merge_name}\n";
}else {
	my $old_name = $nam[0] . $suffix_name;
	my $old_index = $nam[0] . $recal_index_suffix;
	print LOG "$old_name.\n";
	print LOG "$merge_name\n";
	rename $old_name, $merge_name;
	rename $old_index, $merge_index;
	print CHK "${merge_name}\n";
}

my @cal = ($cal, '-i', $intv);
# system $cal, '-i', $intv;
(system @cal) ==0 or die "error: @cal, $?";
print CHK "cal\n";

(system $ann) ==0 or die "error: $ann, $?";
print CHK "ann\n";

my $sum_bam = "1.realn.recal.bam";
my @sum_bam = ($summ, '-s', $sapl_name[0], $sum_bam, $intv);
# system $summ, '-s', $sapl_name[0], $sum_bam, $intv;
(system @sum_bam) ==0 or die "error: @sum_bam, $?";
print CHK "summ\n";

(system $fm_top) ==0 or die "error: $fm_top, $?";
print CHK "top\n";

# my $xls = "1.xls";
my @fm_jxl= ($fm_jxl, "1.xls");
(system @fm_jxl) ==0 or die "error: @fm_jxl, $?";
print CHK "jxl\n";

my $final_name = $sapl_name[0] . "_R12.fq.var.xls";
$final_name = lc $final_name;
rename "1.jxl.xls", $final_name;

if ($cnvcal) {
	print LOG "CNV calling";
	my @quandico = ($quandico, '-n', $sapl_name[0], $cnvcal_sample, $cnvcal_ref, $intv);
	(system @quandico) ==0 or die "error: @quandico, $?";
	print CHK "cnv\n";

}

my $dir_current = ".";
opendir my $dir_cur, $dir_current or die "Error: cannot open ${dir_current}: $!";

foreach (readdir $dir_cur) {

	if ($_ ne '.' and $_ ne '..') {
		rename $_, ".$_";
	}
	
	
}
closedir $dir_cur;

rename ".1.realn.recal.bam", "1.realn.recal.bam";
rename ".1.snp.vcf", "1.snp.vcf";
rename ".$final_name", "$final_name";
rename ".1.dbnsfp.snpeff2gatk.vcf", "1.dbnsfp.snpeff2gatk.vcf";
if ($cnvcal) {
	rename ".1.csv", "${sapl_name[0]}_cnv.csv";
	rename ".1.pdf", "${sapl_name[0]}_cnv.pdf";
	rename ".1_overview.pdf", "${sapl_name[0]}_cnv_overview.pdf";
}


close LOG;

	





















